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Abstract: A method is presented for estimating subaqueous integrated backscatter from the 
CALIOP lidar. The algorithm takes into account specular reflection of laser light, laser 
scattering by wind-generated foam as well as sun glint and solar scattering from the foam 
Analyses show that the estimated subaqueous integrated backscatter is most sensitive to the 
estimate of transmittance used in the atmospheric correction, and is very insensitive to the 
estimate of wind speed used. As a case study, CALIOP data over Tampa Bay were 
compared to MODIS 645 nm remote sensing reflectance, which previously has been shown 
to be nearly linearly related to turbidity. The results indicate good correlation on nearly all 
CALIOP clear-free dates during the period 2006 through 2007, particularly those with 
relatively high atmospheric transmittance. When data are composited over the entire period 
the correlation is reduced but still statistically significant, an indication of variability in the 
biogeochemical composition in the water. Overall, the favorable results show promise for 
the application of satellite lidar integrated backscatter in providing information about 
subsurface backscatter properties, which can be extracted using appropriate models. 

Keywords: lidar; CALIOP; CALIPSO; turbidity; water quality; backscattering, Tampa Bay. 
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The U.S. Environmental Protection Agency [1] ranked siltation and its attendant turbidity as one of 
the most widespread pollutants in the United States. It affects rivers, streams, lakes, reservoirs, 
estuaries, and coastal waters, adversely affecting aquatic habitat, drinking water treatment facilities, 
recreational use, and commercial navigation and fisheries. High levels of sediment in the water column 
reduce spawning grounds and reduce the availability of food supplies for aquatic life. High 
concentrations of sediment may lead to rapid deposition, causing changes in coastal or fluvial 
morphology, which may lead to navigational hazards or increase flood risks for the region. Regions of 
low turbidity due to low suspended sediment concentrations) may experience rapid erosion leading to 
entrenchment and bank failure in a fluvial system or to beach erosion and barrier island collapse in the 
coastal environment. 

Remote sensing of turbidity has a long history in the scientific literature. Descriptions of qualitative 
analyses relating Landsat Multispectral Scanner System (MSS) surface radiance to turbidity and 
suspended sediment concentration appear in the literature as early as 1973 [e.g. 2, 3, 4]. Quantitative 
algorithms began to appear in the 1970s, and ranged in complexity from linear relationships between 
MSS Band 5 radiance and suspended sediment concentration to polynomial fits to the ratio of MSS 
radiances in different bands [see citations in 5]. More rigorous, model-driven algorithms were in 
development by the late 1970s [e.g. 5], with the incorporation of diffuse reflection models developed 
earlier in the decade [e.g. 6, 7, 8]. Curran and Novo [9] summarized early efforts, though several 
significant studies post-date this work [e.g. 10, 11-14]. More recently, good results have been obtained 
in retrieving turbidity using off-nadir Moderate Resolution Imaging Spectroradiometer (MODIS) 
imagery[e.g. 15]. Results indicate a near linear relation between in situ turbidity and red reflectance 
when applied to individual regions. Fonnulation of a global algorithm relation between turbidity and 
reflectance has not been demonstrated. 

Spaceborne lidars have been used in earth observation since the mid-1990s when BALKAN flew on 
the space station MIR [16]. A large fraction of the lidar pulses emitted by these lidars have fallen on 
water. Unlike their airborne counterparts, spaceborne lidars have not had the vertical resolution 
necessary to profile the water column. This paper offers to add value to existing lidar datasets by using 
the depth-integrated lidar backscatter to extract infonnation about turbidity, or scattering within the 
water column. 

The use of a 532 nm laser to study water turbidity was pioneered in the 1980s with theoretical work 
by Gordon [17], and by Phillips et al. [18], Experiments into the use of a ship-mounted laser (532 nm) 
for measuring the optical properties of seawater were carried out by Ivanov et al. [19] in the Soviet 
Union. Empirical studies using the Australian Weapons Research Establishment Laser Airborne Depth 
Sounder (WRELADS) sensor by Phillips et al. [20], and by Billard [21] showed that estimates of both 
backscattering and effective attenuation were possible from depth sounding lidars. The 532 nm 
channel is of particular interest in both turbidity monitoring and in other hydrologic applications, 
because it lies very close to the wavelength achieving maximum penetration into the water column 
[ 22 ]. 

Spaceborne lidars have been used very little for oceanography. Lancaster et al. [23] used the GLAS 
(Geoscience Laser Altimeter System) lidar on the ICESat (Ice, Cloud, and land Elevation Satellite) 
platfonn to detect water surface reflectance and compared that with the predicted theoretical 
reflectance from QuikSCAT (Quick Scatterometer) wind speeds. Hu et al. [24] used the CALIOP 
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(Cloud-Aerosol Lidar with Orthogonal Polarization) lidar on the CALIPSO (Cloud-Aerosol Lidar and 
Infrared Pathfinder Satellite Observations) satellite to measure wind speed, and compared the results to 
AMSR-E (Advanced Microwave Scanning Radiometer-EOS) wind speed measurements. They used 
the depolarization ratio available from the CALIOP instrument to remove the effect of sea foam from 
the surface measurement. 

The study detailed in this paper also uses the CALIOP lidar, which has a 30-meter vertical 
resolution near the surface. CALIOP has a 90 m spot diameter, with spots spaced 333 m along track. 
CALIPSO flies in NASA's A-Train constellation, following Aqua by 73±43 seconds. [25] As an active 
sensor, it is capable of making night-time measurements, and the ability to separate atmospheric 
returns from surface returns make local atmospheric correction possible. 

The overall approach is to estimate subaqueous backscatter as the residual term within an 
expression of total depth-integrated attenuated backscatter. Specifically, total, depth-integrated 
attenuated backscatter is formulated as the sum of all surface and subsurface effects. With total 
backscatter obtained from CALIOP, and other surface components modeled, the subsurface 
backscatter can be retrieved. Sensitivity analysis on the accuracy of the CALIOP data and on 
individual model components is used to delineate the limits potential accuracy of the subsurface 
backscatter retrieval. The theoretical sensitivity analysis is then followed by a test case over the 
Tampa Bay region of Florida. 

2. Model basis 

The depth-integrated attenuated backscatter (at a wavelength of A, in nm) received by a satellite 
from the surface bins can be represented as an attenuated sum of scattering returning from components 
of the target, or 

Yx = T a (y7 + yI + y! + Yx +Yx) ( 1) 

T) is the atmospheric transmittance (including aerosol effects, and along the lidar look direction), the 
variables Yx’Yx’YL Y;/\ and Yx represent (respectively) the integrated backscatter due to: specular 
reflection of the laser from the water surface, laser light scattered by foam on the water surface, the 
specular reflection of sunlight from the water surface, sunlight scattered by foam on the water surface, 
particles in the underwater environment, y'x , y[ - Y~h and Yx ^ are all primarily controlled by the wind 
and wave field. Yx ' s controlled by the number and properties of particulates present in the water 
column. Other possible components not listed in (1) include bottom scattering effects, where the water 
is shallow and clear enough that photons can reach the sea floor and scatter back to the sensor; 
scattering from surface slicks, where they exist; and scattering by spray in regions of high winds. 
These factors are beyond the scope of this study. Details of each component are described below. 

The goal of the modeling presented here is to allow the estimation of Y 53 2 from satellite 
measurements of y 53 2 , and y 1064 , an d the modeling of the other components. Surface wind speed and 
visibility estimates are used as inputs into the models. 


2.1. Atmospheric Correction 
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A principal factor affecting the radiometric surface measurement is the atmospheric transmittance, 
which can be further subdivided into molecular, aerosol, and cloud transmittance. One advantage of 
using lidar technology compared with a passive sensor is the ability to separate the scattering due to 
the atmospheric effects from the scattering at the surface based on the time of the photons' arrival at 
the satellite. The CALIPSO release products include an aerosol layer optical depth and a cloud layer 
optical depth product. These data, combined with temperature and pressure data from NASA’s Global 
Modeling and Assimilation Office, and included in the CALIPSO release product, were used in the 
current work to develop a LOWTRAN atmospheric model. 

In the current study, optical depth of the cloud layers and aerosol layers were accounted for using 
the CALIOP data. However, optical depths of the clear air interstices, not available from the CALIOP 
archive, were evaluated using the LOWTRAN model, resulting in a continuous profile of optical depth 
throughout the atmospheric column. The optical depths for the layers were summed (taking into 
account that the aerosol and cloud layers sometimes overlap), and a resulting total atmospheric optical 
depth and transmittance were obtained. To calculate the clear-air optical depths, default settings were 
used for the LOWTRAN model, except for the input of relative humidity, temperature, pressure and 
ozone concentration data extracted from the Level 1 and Level 2 CALIOP products, and surface 
visibility was interpolated from historical airport records from around Tampa Bay. A sensitivity 
analysis showed that assuming a constant surface visibility of 23 km produced errors in the subaqueous 
integrated backscatter of less than 5%. 

2.2. Specular surface reflection 


The specular reflection of the lidar from the water surface is by far the largest contribution to the 
total integrated attenuated backscatter. The specular reflection decreases as wave steepness increases, 
and is extremely large in very calm waters. 

Hu et al. [24], citing Platt [26], Menzies et a/.[27], and Tratt et a/.[28], modeled the specular surface 
reflection, where px is the Fresnel specular reflection coefficient (p 532 * 0.0209, p w 64 « 0.0199), o 2 
is the wave slope variance (see Equation 3), and 6 is the zenith angle of the sensor: 


7 " = 4 no 2 cos*9 eXP 


-tan 2 6 1 


Here, the small angle approximation is not made, as it was made by Hu et al., because it can 
produce significant errors at low wind speed even for small angles such as the 0.3° pointing angle for 
CALIOP. Hu et a/.’s composite model for the wave slope variance as a function of wind speed is 
employed: 

f U < 7 0.0146VZ7 

o- 2 = < 7 < U < 13.3 0.003 + 0.00512 1/ (3) 

( U> 13.3 0.138 log 10 U - 0.084 


Although this model provides a reasonable guess at the backscattering when given a known wind 
speed, this model shows strong sensitivity to the accuracy of that known wind speed. To overcome this 
sensitivity, one can exploit the fact that the only wavelength dependent term in Equation 2 is p A , and 
the fact that infrared light does not penetrate the ocean surface; i.e. 
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7 io 64 is zero. We can then estimate the specular surface reflection in the green, using the total 
integrated attenuated backscatter in the infrared (71064) : 


w P 532 w P 532 (Y 1064 

7532 = 7 71064 = 


Pl064 


Pi 064 ^1064 


7l064 


7l064 " 


7io / 6 4 ) 


(4) 


7/0 64 , 7(064, and 7 1( , 64 arc modeled as described in the subsequent sections. The resulting model is 
only weakly dependent on the accuracy of the wind speed estimate, especially for low wind speeds and 
for night measurements, when yl 064 and Ywm are zero. 

An added benefit of this infrared-based approach is that deviations from the fully-developed sea 
models of Cox and Munk [29] or Breon and Henriot [30] are accounted for implicitly, without the need 
for estimating fetch, wind direction, or bottom topography. This is particularly important in coastal 
zones, where sheltered waters can have very short fetches, and shallows can cause rapid and localized 
steepening of waves. 


Figure 1. Data of Frouin et al. [31], which suggests an exponential function of wavelength. 
The correlation is not significant, with only four data points, but it suggests that 
exponential decay is a starting point upon which we can construct a model. 



Wavelength (nm) 


2. 3. Lambertian scattering from foam 

At relatively high wind speeds, the scattering of the lidar beam from whitecaps and foam streaks on 
the water surface can be a significant factor, particularly as the wave surface return becomes darker as 
the waves steepen. At wind speeds higher than about 10 m/s, the scattering from foam begins to fall 
within the range expected for subsurface scattering in low turbidity natural waters. As noted above, Hu 
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et al. [24] used the depolarization ratio from the CALIOP lidar to model the foam contribution to the 
backscattering. However, depolarization may occur during scattering from foam or from hydrosols, so 
this type of model may eliminate the potential to retrieve information from the subsurface. The foam is 
therefore modeled as a Lambertian scattering process, from foam covering a fractional area described 
by Callaghan and White [32]: 

( U < 3.70 0 

W = 3.70 < U < 10.1874 3.18 x 10~ S (U - 3.70) 3 ( 5 ) 

l U > 10.1874 4.82 x 10 ~ 6 ([/ + 1.98) 3 

Moore et al. [33] modeled the reflectance of foam as a function of wind speed. In this model the 
reflectance of the foam is expressed as an “additional” reflectance, representing the increased 
reflectance of the ocean surface due to the foam. In the visible wavelengths, including 532 nm and 
expressed at the upper limit of 670 nm, they express this additional reflectance as: 

RSAR< 670 = 3.14 x l(T 6 t/ 2 - 55 . (6) 

For the near infrared at 860 nm they express the additional reflectance, RSARxeo, as 

RSAR 860 = 0.22 (1 - exp [—4.2 RSAR 670 ]). (7) 

Figure 2. Plot showing fit of Equations 9a (black) and 9b (red) to points created by fitting 
an exponential decay for each of several constant wind speeds. The fits shown create 
functions of wind speed that can be used in Equation 8. 
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The data of Frouin et al. [31] (See Figure 1.) suggests an exponential decrease in the reflectance of 
foam with wavelength from laboratory tests. 
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RSAR a>670 = A[U]exp [~k[U]A]. (8) 

Figure 3. Model of backscatter due to Lambertian scattering from foam. Three 
wavelengths are shown, the black points show the model of Moore, et al. [33], for visible 
wavelengths, up to 670 nm (Equation 6). This model is used for the 532 nm laser light. The 
blue points show the model of Moore, et al. [33] , for 860 nm shortwave infrared (Equation 
7). The red curve shows an extrapolated model for 1064 nm, created (as described in the 
text) as an exponential decay with wavelength perfectly fitting the two shorter-wavelength 
models. (Equation 8, where 2 = 1064 nm.) 



♦ Visible, up to 670 nm ©860nm ♦1064 nm 


At a given wind speed, the coefficient A and the decay constant k can be estimated by fitting to 
equations 6 and 7. Repeating this process over a range of wind speeds, and using least squares 
polynomial regression, we can estimate A[U] and k[U] as follows. 

A[U] = 1.53 x 10“ 4 + U (-1.17 x 10“ 4 + {/( 2.57 x 10“ 5 + [/(- 2.27 x 1(T 7 + U 1.74 x 10“ 8 ))) (9 a ) 

k[U] = 4.16 x 10 5 + y (-3.02 x 10 2 + U (9.86 x 10 1 + U{ 5.30 + [/(- 2.68 x 1(T 2 )))J (%) 

Using these results, and multiplying by cos ® / n to convert between reflectance and Lambertian 
scattering, one can write the expression for the scattering due to the foam. Because the model 
presented above calculates the additional reflectance due to the foam, we also add in an estimate of the 
sea surface scattering without foam, for which we use Equation 2. It is possible to use Equation 10b for 
the 532 nm lidar, and the same results will be obtained. Equation 10a can be used when efficiency is a 


concern. 
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yL = w I 


P 532 


4 n er 2 cos 4 0 


Pl064 


-exp 


-tan 2 6 


2er 2 


+ 3.14 x 10 ~ 6 U 2 - 55 


cos 6f 


n 


^ = f 4,a^ eXP 


— tan 2 0 


2er 2 


+ 21 [17] exp [- 1064/c [t/]] ■ 


cos 0 N 


7T 


(10a) 

(10b) 


2.4. Sun glint 


At low solar zenith angles, and particularly for lidar sensors which may use low-power micropulse 
sensing, the presence of specular reflection of solar radiation, or sun glint, may affect the 
measurement. To estimate the sun glint, the following model, after Kay et al. [34] is employed. 


Yl = (1 - wo 


E^p\(^, 3i]p\ z x> Zy\ ^,1-cos e s sec o v 

4 71 cos 4 /? x 


( 11 ) 


where W is the fractional area covered by foam, Ex is the extra-atmospheric solar irradiance, p(co,X) 
is the Fresnel reflection coefficient of the water surface for the given solar/sensor geometry, and 
p(z x ,z v ) is the probability of a given surface element having the required surface slope to specularly 
reflect light from the sun into the detector, 6 is a zenith angle, and, in subscript, s indicates a solar 
angle and v indicates a sensor angle, to and [3 are angular measures accounting for the geometry of the 
sun and sensor: (where <p is an azimuth) 


cos 2 a) — ^(sin 6 S sin 6 V cos(<p s — q) v ) + cos 6 S cos 9 V + 1) 


cos /? — 


cos 9 S + cos 9 V 
2 cosa) 


(12a) 

(12b) 


Kay et al. [34] provide the Cox and Munk [29] expression for p(z x ,z v ) as well as a more recent 
version by Breon and Henriot [30] which presents coefficients with reduced uncertainty. In this study, 
the version of Breon and Henriot is used. 

In the absence of foam, sun glint, while dependent on wind speed through the probability density 
function for wave slopes, is relatively insensitive to wind speeds over the normal range, changing only 
by 0.3% between 0 and 30 m/s. The error in the scattering estimate caused by sun glint is also very low 
compared with the intensity of the lidar backscattering, because only approximately 0.006% of the 
sun’s energy falls within the bandwidth of the lidar detector. In the presence of foam, the sun glint falls 
off as the wind speed increases because less glint is possible from foam-covered slopes. This model 
assumes that foam is distributed over the surface in a pattern uncorrelated with water surface slope. 
Even in the case where this assumption is invalid, the estimated sun glint for 3.7 m/s represents an 
upper bound on the sun glint for wind speeds greater than this value, as it represents a foamless case. 
For wind speeds less than 3.7 m/s, this assumption is not necessary, as foam is assumed to be absent. 

The sun glint model is relatively insensitive to uncertainty in the parameters. Using the published 
[30] uncertainties in the parameters for the wave slope probability, and 10% uncertainties in the 
parameters of the foam model as above, the model exhibits, at most, a 3% error in the predicted sun 
glint. 
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2.5. Sunlit foam 


The contribution due to sunlit foam is modeled making the assumption that foam is a Lambertian 
scatterer. Using a method similar to Moore et al. [33], the scattering is modeled as an additional 
scattering tenn dependent upon the foam reflectance. 


Yf 


a- 


(13) 


The first term in this equation represents the surface return in the absence of foam, while the second 
term represents the additional scattering due to the foam. W is the area covered by foam (Equation 5), 
Yx is the sunglint scattering (Equation 11), RSAR, is the additional foam reflectance (Equation 6 or 8, 
depending on wavelength), 6 is a zenith angle, and, in subscript, 5 indicates a solar angle and v 
indicates a sensor angle, Ex is the extra-atmospheric solar irradiance, Q is the solid angle of the lidar 
detector, A is the area of the ground spot, / is the laser energy, t p is the pulse length, and T, is the 
atmospheric transmittance (along the lidar look direction). See Table 1 for the lidar-specific parameters 
applicable to CALIOP. The only wind-speed dependence for this term is the areal coverage of the 
foam. For all wind speeds less than 30 m/s, the foam contribution is an order of magnitude less than 
the sun glint. 

The sensitivity to model parameters in the sunlit-foam model is dependent on wind speed. The 
parameters of the sunlit foam model are the same as those of the lidar foam model, and relate to the 
coverage and reflectance of the foam. However, the sunlit-foam model is less affected by the model 
parameters than the lidar foam model. At low-to-moderate wind speeds, the model is fairly insensitive 
to the model parameters, with uncertainties less than a factor of two for wind speeds less than about 
22 m/s. 


Table 1. CALIOP lidar-specific parameters for Equation 13. 


Parameter 

Value 

Note 

O 

1.26 x 10" 13 sr 


A 

3.85 x 10 3 m 2 


1 

1.1 x 10 1 J 


tp 

2.0 x 10" 8 s 


o v 

5.236 x 10' 3 (0.3°) 

(June 2006 - Nov 2007) 

5.236 x 10' 2 (3.0°) 

(Nov 2007 - ) 


2. 6. Subaqueous scattering 

Solving for the subaqueous scattering (y^ 2 ) using Equation 1, and combining with Equation 4, the 
full model can be expressed as: 

u 7532 ( P 532 (Y 1064 / c Sf \ , f , S , 

Y532 = r I y 71064 - KiW - 7 i 064 ) + 7i +Ya+ 7/ j (14) 

‘ 532 1064 '1 1064 ' ' V ' 

Combining Equation 14 with Equations 6, 8, 10, 11, and 13, we can expand this into the following 
fonn: 
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V 532 — ' 


VS32 


1 532 


t 


P 532 / 71064 


Pl064 \ ^1064 


- W 


P 1064 


4 7i a 2 cos 4 9 


exp 


-tan z 6 


2(7 2 


+ A[U ] exp [ — 1064 /c [LA]] - 


cos 6 




- (1-M0 


E 1064 p[(*>, 1064]p[. 


z x> z y\ rril—cos 0 S sec 0 V 


4 71 cos 4 /? 


'1064 


yy I Exp\u>, 1064 \p\z x ,z y r ^i- cos g sSec g v 


4 7i cos 4 ^5 
+ W(A[U]ex p [-1064 k[U]]) 

PS32 


cos 8 S E 1064 UA cos 9 v i - cos e sec e, 


+ W 


4 7T <J 2 COS 4 # 


exp 


cos 0, 
—tan 2 # 


/ 


2cr 2 


+ 3.14 X 10 1 / 


6 jt2.55 


cos 9 


I I c-t T/l ,^ E^ 32 p[<n, 532]p[z x , Zy ^i-cosg.jsecgp 
+ |U ; 4ticos 4 /? 532 

I | I/I t (E532P ^ 32 ]p[z x ,Zy] rr i- cos g sSec g v 

+ l 1 4 71 cos 4 ^5 532 


rr COS 0 ^ EconflA COS 0 V -1 a car ft 

+ W ( 3.14 X 10 “ 6 t / 2 5S ) - 71 _ cos fleece,, 

COS#,, I 532 


\ 

) 

/ 


7T 


(15) 


For nighttime overpasses, sun glint and sunlit foam (y|, y7 ) in both wavelengths are zero, so this 
simplifies to Equation 16. 

„ K532 

K532 — 



- iy 


+ iy 


+ i4[t/]exp [— 1064fc[t/]] ■ 


cos # 


7T 


+ 3.14 X 10 ~ 6 U 2 


cos #^ 


(16) 


2. 7 . Comparison of backscatter terms 


Figure 4 shows model results for the magnitude of the components of the total integrated attenuated 
backscatter (excepting subsurface scattering), to give a sense of the relative importance of each term, 
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assuming a solar zenith angle of 30°. The specular wave scattering shown is from Equation 2, in 
contrast to the method outlined in Section 2.2, as this allows a theoretical dependence on wind speed, 
for comparison to the other terms. In the case study, subsurface scattering ranged from 0.00 (or below) 
to 0.13 sr' 1 . 

Figure 4. Theoretical CALIOP scattering off water surface as a function of wind speed . 

Wave scattering is modeled using Equation 2. Foam scattering and sun-on- foam are 
modeled by Equation 10. Sun glint is modeled as Equation 11. Graphs illustrate the relative 
importance of the different components to the total integrated scattering. The solar zenith 
angle is assumed to be 30° when computing the sun glint and sunlit foam components. 
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2.8. Uncertainties 

Uncertainty propagation analysis allows the evaluation of the sensitivity of the model to input error. 
There are five major inputs into the model: First, an estimate of wind speed is required to estimate the 
foam component. Second, transmittances are needed for both the green and the infrared channels. 
Third, the total integrated attenuated backscatter from the lidar measurements is needed, in both green 
and infrared; uncertainty is introduced both during the measurement process and also during any 
numerical integration with depth. 


532 nm Foam 

— • - 1 064 nm Sunglint 
532 nm Total 



2.8.1. Wind Speed Estimate 
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The use of Equation 4, instead of Equation 2 for the specular surface return results in a modeled 
subaqueous integrated backscatter that is relatively insensitive to the wind speed estimate. 

Figure 5. Error in subaqueous integrated backscatter due to 1 m/s overestimate in wind 
speed estimate. Because the wind speed estimate is only used to estimate foam cover and 
reflectance, the error at low wind speeds is negligible. 



2.8.2. Transmittance Estimate 


Two transmittance estimates are used in the atmospheric correction of the lidar data; at 532 nm and 
at 1064 nm. The accuracies of these estimates are important considerations in determining the accuracy 
of the subaqueous integrated backscatter estimate. From Equation 15, it can be shown that: 


dy 532 _ 2 ^ 532 

^532 T 5 3 32 


(17a) 


532 _ 2 PS32 V 1064 
d7l064 Pl064 Ti 06 4 


(17b) 


Equation 17 shows that the sensitivity to the transmittance increases linearly with the total 
backscatter. Figure 6 shows the dependence upon both the total integrated backscatter and the 
transmittance. It is clear from both Equation 17 and from Figure 6 that the accuracy of the 
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transmittance estimate plays an important role in determining the uncertainty of the final result, and 
that low transmittance data may provide largely meaningless results. 

Because of the opposing signs of (17a) and (17b), positively correlated errors in T 532 and Tn )64 will 
tend to cancel. This will tend to limit the error introduced by the surface visibility parameter to the model. 

Figure 6. Error in subaqueous integrated backscatter due to 0.01 overestimate in the 532 
mn transmittance estimate. The error due to an overestimate in the 1064 mn transmittance 
is proportional, with a constant of proportionality oi — p WM / ~ —0.9494. The error 
is dependent upon the total integrated backscatter received; the two curves represent low 
(0.01 sr , blue) and high (0.1 sr' 1 , red) integrated backscatters. 



2.8.3. Integrated Backscatter 

The attenuated backscatter measured by the lidar sensor is numerically integrated to obtain the total 
integrated attenuated backscatter. There is uncertainty associated both with the original attenuated 
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backscatter measurements and with the numerical integration scheme used. The sensitivity of the 
model to these uncertainties can be summarized by the following: 


5*32 _ 
^7532 


(18a) 


dy 532 _ P 532 rr _2 
— — ~ ' 106 
°Y 1064 Pl064 


(18b) 


Like the sensitivity to the transmittance, the sensitivity to the total integrated attenuated backscatter 
is dependent on the transmittance. The error in the subaqueous integrated backscatter is always greater 
than the error in the total integrated attenuated backscatter, and is greater by more than a factor of four 
for transmittances less than 0.5. This sensitivity makes the accuracy of the total integrated attenuated 
backscatter input at least as important as the accuracy of the transmittance estimate in determining the 
uncertainty in the final result. 


Figure 7. Error in the subaqueous integrated backscatter for a 0.001 sr-1 error in the total 
integrated attenuated backscatter. This level of uncertainty was observed in the CALIOP 
data when comparing different numerical integration schemes. The errors associated with 
uncertainties in both bands are shown. 
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As for the transmittances, the opposing signs (and roughly equal magnitude; p 532 /Pio64 ~ 1-053) of 
(18a) and (18b) will cause the effects of positively correlated errors in y 532 and yio 64 to cancel. This 
should help to ameliorate the effect of the choice of numerical integration scheme, which might 
introduce systematic errors. 
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3. Case Study: Tampa Bay, Florida, USA 

3. 1. Site Description 

Tampa Bay is located on the west coast of the Florida Peninsula between 27.3-28. 1°N and 82.4- 
82.8°W (Figure 8). Fresh water inflow to the bay comes from a watershed of about 596000 km“ [35], 
through four major rivers - all of them controlled: the Hillsborough River, the Alalia River, the Little 
Manatee River, and the Manatee River. Chen, et al. [15] found that the Alalia and the Hillsborough 
dominate the delivery of dissolved organic carbon (CDOM) to the bay, and also found that CDOM is 
the primary absorber (at blue wavelengths) in the estuary. Lewis and Whitman [36] identified seven 
major segments of Tampa Bay: Old Tampa Bay, Hillsborough Bay, Middle Tampa Bay, Lower Tampa 
Bay, Boca Ciega Bay, Terra Ceia Bay, and Manatee River. In general, sediment size increases as one 
progresses from Hillsborough Bay through the Lower Tampa Bay. Old Tampa Bay is heterogeneous, 
though generally it coarsens toward the outlet. Terra Ceia, and the Manatee River estuary are 
characterized by medium to fine sands. Boca Ciega Bay contains fine sands to muds, generally 
coarsening to the south. [37] 

The Environmental Protection Commission of Hillsborough County (EPCHC) samples turbidity 
and other water quality parameters at 56 fixed stations in Tampa Bay once per month, on a rotating 
schedule. 

3.2. Data 

Data from the CALIOP lidar were combined with reflectance data from the MODIS radiometer. 
Surface wind speed and visibility data, from airports and buoys, along with atmospheric data provided 
with the CALIOP data were used to construct a LOWTRAN model. 

3.2.1. CALIOP data 

CALIPSO, the satellite upon which the CALIOP lidar resides, was launched on 28 April 2006. 
Seven overpasses were selected between that date and the late November, 2007 pitch maneuver that 
changed the zenith angle of CALIOP from 0.3° to 3.0°. The seven overpasses selected for evaluation 
were chosen based on absence of cloud cover. The descending pass of CALIOP is the only pass that 
traverses the bay, which forced the use of solely nighttime scenes. Version 3 CALIOP data were used 
except where noted. 

Four data products were used in this analysis. The Level 1 product provides backscatter data as a 
function of altitude. For altitudes below sea-level, the altitudes were corrected for the speed of light in 
seawater compared to air, and the resulting function of altitude was numerically integrated using a 
trapezoidal scheme over the interval from the surface bin through five bins below the surface. Low- 
pass filtering inherent in the CALIOP electronics causes the analog-to-digital conversion to smear the 
surface return into subsequent bins, so this integration is necessary, although these five bins represent a 
depth of water (from which it is unreasonable to expect real returns. [24] 

The Level 2 Cloud and Aerosol Layer products were used to provide the optical depths, and the 
altitudes of the tops and bottoms, of the atmospheric strata influenced by these phenomena. Shots that 
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were flagged as containing opaque layers or other quality problems were rejected from the analysis. 
The Level 2 Atmospheric Profile (Version 2 beta) product was used only to provide temperature, 
pressure, ozone, and molecular profile data, all provided by the NASA Global Modeling and 
Assimilation Office (GMAO), which was then used as input into a LOWTRAN model. 

3.2.2. MODIS data 

Terra-MODIS Level 2 (MOD09QK) surface reflectance (nonnalized radiance) products were 
retrieved from the LPDAAC corresponding to the afternoon before each CALIPSO overpass. The 
Terra overpass was a compromise of the closest overpass in time and the closest to nadir observation. 
After geocoding, the data located over Tampa Bay were extracted. The nonnalized radiance values (p ) 
were converted to remote-sensing reflectances (Rrs), assuming that the illumination can be 
approximated as isotropic: 

Rrs = ~ (19) 

n 


3.2.3. Supplementary data 

Historical meteorological data were obtained from seven airports and four buoys situated in and 
around the bay. The airports were: Tampa International Airport (TP A) [27° 58’ 32”N, 82° 32’ 0”W], Peter 
O. Knight Airport (TPF) [27° 54’ 56”N, 82° 26’ 57”W], MacDill Air Force Base (MCF) [27° 50’ 57”N, 82° 31’ 
15”W], Tampa Executive Airport (VDF) [28° 00’ 50”N, 82° 20’ 43”W], St. Petersburg/Clearwater 
International Airport (PIE) [27° 54’ 36”N, 82° 41’ 15”W], Albert Whitted Airport (SPG) [27° 45’ 54”N, 82° 
37’ 36”W], Sarasota/Bradenton International Airport (SRQ) [27° 23’ 44”N, 82° 33’ 15”W]. These airports 
record hourly measurements of wind speed and visibility. The several Coastal-Marine Automated 
Network (CMAN) buoys that provide six- to twenty-minute wind speed data are Egmont Key 
(EGKF1) [27° 36’ 3.6”N, 82° 45’ 36”W], Anna Maria Island (ANMF1) [27° 33’ 0”N, 82° 45’ 0”W], Port 
Manatee (PMAF1) [27° 38’ 13”N, 82° 33’ 47”W], and Clearwater Beach (CWBF1) [27° 58’ 36”N, 82° 49’ 
54”W], which provided data for this study. The locations of these ground stations within the 
topobathymetric setting of the estuary are shown in Figure 8. 

Depths were acquired from the NOAA National Geophysical Data Center (NGDC) Coastal Relief 
products, combined with gridded airborne lidar depth soundings from the Experimental Advanced 
Airborne Research Lidar (EAARL) in some shallow waters; these data were also acquired from the 
NGDC. The EAARL data are limited to depths shallower than three meters, and are not spatially 
extensive. These depth data were gridded together at approximately 3.67 m posting, as this was the 
posting of the processed lidar dataset, to obtain seamless bathymetry of the bay. In most cases, 
however, the resolution of the depth data is significantly less than this posting. 

Although Tampa Bay has a strong history of in situ turbidity research, there are no data both 
coincident and contemporaneous with the lidar coverage in the period of study. Although the small (90 
m) lidar spot size is an advantage when considering spatial heterogeneity, it becomes a liability when 
attempting to find fortuitous matches with established study sites. Over the period between June 2006 
and October 2007, there were no study sites within 1 km of any overpass that were occupied within 
two days of the overpass. 
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Figure 8. Topobathymetry of Tampa Bay. Bathymetry in purple is derived from an 
airborne lidar survey in May 2007 [38]; topography and bathymetry in blue is a composite 
of NOAA GEODAS and the Global Multi -Resolution Topography product [39]. Also 
shown are the seven airports and four CMAN buoys from which wind speed and/or 
visibility data were obtained. The map is shown in a latitude-longitude pseudoprojection, 
with northwest comer [28°04’26.546304”N, 82°52’08.082264”W], and southeast comer 
[27°19’58.28718”N, 82°19’58.004148”W]. 
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3.3. Analysis 

The data were processed through the algorithm described above (though the sun glint and sun-on- 
foam sections of the model are forced to be zero when the sun is below the horizon). Wind speeds 
were estimated by interpolation over space and time from records collected at the airports and buoys 
shown in Figure 1 . Surface visibility, a parameter in the clear-air LOWTRAN model, was estimated by 
similarly interpolating, though only the airport records were available for visibility. A sensitivity 
analysis was carried out on the effect of these inputs, and it was found that changing the wind speed by 
90% and the visibility by 50% had a nearly undetectable result on the resulting subaqueous integrated 
backscatter. This is due to the wind speed estimate affecting only the foam reflectance, which is small, 
and the visibility affecting only the clear-air optical depth, which represents only a very small fraction 
of the total atmospheric optical depth. It is likely that at very high wind speeds or very low visibilities, 
beyond the range of those seen in this study, these sensitivities would increase. 

Fourteen points were identified as outliers, using Peirce’s criterion, as outlined in the work of Ross. 
[40] All of these points had estimated y u values above 1.0 sr' 1 , more than six times the largest non- 
outlier. All of these points were acquired on 07 May 2007, a date which had relatively low atmospheric 
transmittance. 

Although there are no coincident in situ measurements of backscatter or turbidity, Tampa Bay has 
been well studied as a venue for the use of MODIS 645 nm radiance and reflectance data for the 
estimation of turbidity. MODIS is a passive radiometer, which, in regions far from the specular solar 
disk, images mainly subsurface scattering and relatively low levels of foam and sun glitter from waves. 

In order to evaluate the ability of the CALIOP to detect subsurface scattering, we compare the 
CALIOP data with MODIS Channel 1 data. There are important limitations to this comparison. 
Perhaps the most important is that the MODIS acquisitions require daylight, which imposes a 
minimum time separation between the MODIS acquisition and the CALIOP acquisition; for the seven 
days used, the Terra overpass from the afternoon before the CALIOP measurement was used for 
comparison. During daylight hours, the A-Train configuration of the Aqua and CALIPSO satellites 
could be exploited to derive nearly simultaneous measurements. Another limitation to this comparison 
is the unknown penetration depths for MODIS and CALIOP. Chen, et al. [15] and Moreno-Madrinan, 
et al. [41] showed that the MODIS Channel 1 (645 nm) remote sensing reflectance (Rrs) is well 
correlated with turbidity in Tampa Bay, for depths greater than about 2.5 m (Chen, et al.: 2.8 m; 
Moreno-Madrinan, et al.: 2.4 m). This suggests that the penetration depth for MODIS is in this range, 
and that for shallower waters, the bottom scattering becomes an important component to the total 
return. No previous work has investigated the penetration depth of CALIOP; this depth must depend 
upon the strength of the light, the sensitivity of the detector and the strength of the backscattering from 
the target, so it is not trivial to determine. An investigation of the correlation of CALIOP backscatter 
with MODIS reflectance showed no threshold depth in the relationship similar to those observed 
previously between MODIS reflectance and turbidity. 

It would be preferable to compare the CALIOP data with field observations, as was done with the 
MODIS data by Chen and Moreno-Madrinan. However, the limited spatial extent makes this extremely 
difficult, as only up to five of the EPCHC field sites are within one kilometer of any CALIOP track. 
The time elapsed between measurements typically varies between two days and nine days. Without a 
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planned simultaneous observation campaign, such field validation is not feasible with a single-beam 
lidar. 

Following Moreno-Madrinan, et al. [41], the MOD09QK (quarter-kilometer) land surface 
reflectance product was used for the evaluation. For each CALIOP observation over water, the closest 
MOD09 grid cell was chosen which was entirely water (>99.9%). The NOAA Global, Self-consistent, 
Hierarchical, High-resolution Shoreline Database (GSHHS) [42] was used to identify pixels (both 
CALIOP and MODIS) contaminated by land. The identification of land-contaminated pixels is limited 
by the uncertainties in geolocation of the two instruments; the MODIS uncertainty is increased because 
gridded data were used in this study, and precise geolocation information is not carried into the grid 
format. 

Table 2. Paired CALIOP and MODIS observations of Tampa Bay from June, 2006 
through November, 2007. Excluded from these analyses are fourteen points from 07 May 
2007 that were identified by Peirce’s criterion as outliers. (All of these had y u >1.0 sr' 1 ). 


Date/UTC 
Time of 
CALIOP 
Overpass 

Number 
of Cloud- 
Free 
Shots 

Mean 

Atmospheric 

Transmittance 

(Cloud-Free) 

Date of Terra 
MOD09 
(Daily) 
Product 

Number of Cloud- 
Free Pairs of 
CALIOP/MODIS 

Correlation 
Coefficient between 
CALIOP and MODIS 
(95% c.i.) 

p-value 

2006-08-08 

07:21:50 

144 

0.66 

2006-08-07 

92 

0.26 (0.12; 0.42) 

1.8 x 10~ 7 

2006-09-25 

07:17:41 

148 

0.65 

2006-09-24 

136 

0.23 (0.11; 0.36) 

3.4 x 10' 9 

2007-05-07 

07:23:10 

140 

0.46 

2007-05-06 

70 

0.03 (0.002; 0.16) 

1.2 x 10 1 

2007-05-23 

07:24:13 

160 

0.61 

2007-05-22 

63 

0.43 (0.24; 0.60) 

5.7 x 10' 9 

2007-07-10 

07:24:00 

147 

0.66 

2007-07-09 

113 

0.10(0.02; 0.22) 

8.2 x 10 4 

2007-09-28 

07:16:58 

153 

0.70 

2007-09-27 

53 

0.14(0.01; 0.34) 

6.6 x 10' 3 

2007-10-14 

07:14:53 

143 

0.68 

2007-10-13 

133 

0.46 (0.33; 0.58) 

2.7 x 10 19 

Aggregated 

Data 

1035 

0.64 

Aggregated 

Data 

660 

0.11 (0.07; 0.16) 

5.8 x 10 19 


The seven paired CALIOP and MODIS datasets are shown in Table 2. Also shown are the 
correlation coefficients and significance p-values for each pair. The final row in the table shows the 
results of the analysis for all of the data combined. 
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Figure 9. MODIS- (MOD09-) derived 645 nm Remote Sensing Reflectance (Rrs) (shades 
of blue and cyan) for 13 October 2007 and the corresponding 14 October 2007 CALIOP 
shots (red dots). Land was masked using the NOAA Global Self-consistent, Hierarchical, 
High-resolution Shoreline Database (GSHHS). In this figure, it has been filled with the 
Global Multi-Resolution Topography Product. The very high values of Rrs along the 
shorelines are probably due to bottom reflection. The CALIOP footprints are shown larger 
than actual size for visibility. 
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Figure 10. Transect across the Tampa Bay Region (see Figure 2 for location) following the 
14 October 2007 CALIOP track, a) CALIOP integrated subsurface backscatter (y u ). b) 
MODIS remote sensing reflectance (Channel 1: 645 nm). The MODIS data presented 
represents the closest MODIS grid cell (to a CALIOP pixel) which is free of land, c) 
Transect of elevation/depth along the CALIOP track, relative to the NAVD datum. 
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3.3. Discussion 

As an example, the MODIS Rrs of 13 October 2007, calculated from the MOD09 Land Surface 
Reflectance product, is shown in Figure 4, where the land (from the GSHHS database) is masked out 
and filled with topography. Superimposed on the MODIS reflectances are the CALIOP shot locations 
from the overpass of the early morning (local time) of 14 October 2007. Figure 10a shows the results 
of the subsurface backscatter algorithm along the CALIOP swath. The MODIS Rrs along the swath is 
shown in Figure 10b. This corresponds to the data in Figure 2, along the CALIOP track. Figure 10c 
shows the depth profile along the track, as a reference for the other two plots. Four distinct zones have 
been identified on this plot; from south to north: Gulf of Mexico is the area southwest of Tampa Bay, 
offshore from Anna Maria Island. It is characterized by low reflectance in the MODIS 645 nm 
channel. Lower Bay is the stretch of water between Anna Maria Island and Mullet Key, which is the 
mouth of the bay, spanning the outlet of the navigation channel. It is characterized by moderate to low 
reflectance in the MODIS 645 nm channel. Boca Ciega Bay is a region of shoals and grass beds 
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between Mullet Key and the mainland, southwest of St. Petersburg. Boca Ciega Bay is characterized 
by high reflectance and high scatter in the 645 nm channel, and is characterized by extremely shallow 
waters (< 4 m). Old Tampa Bay separates St. Petersburg from Tampa, and is characterized by shallow 
waters ( < 5 m) and can have a wide range of reflectance values, though when it has high reflectance, 
the values are more self-consistent than Boca Ciega Bay. 

Figure 1 1 is a scatterplot of the collected CALIOP y u and the MODIS 645 nm Rrs. The correlation 
is not particularly good, (Pearson’s correlation coefficient is 0.11 (0.07; 0.16).) but it shows a 
remarkably high level of significance. It is clear that in the mean, the CALIOP and the MODIS data 
are responding to similar phenomena.. 


Figure 11 . Scatter plot of CALIOP 532 nm integrated subsurface backscatter and MODIS 
645 nm remote sensing reflectance (as derived from the MOD09 Land Surface Reflectance 
Product). The Pearson’s correlation coefficient is 0.11 (0.07; 0.16). The exceptionally low 
points {y u < -0.05 sr' 1 ) are all from 07 May 2007, from points of low transmittance. 
Fourteen outliers were identified with Peirce’s criterion and are not shown. (All of these 
had y u >1.0 sr' 1 .) 
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4. Conclusions 

A method has been presented for the estimation of subaqueous integrated backscatter from the 
CALIOP lidar. External inputs to the algorithm are limited to wind speed and surface visibility, and to 
neither of these inputs does it exhibit high sensitivity. The algorithm takes into account specular 
reflection of laser light at the water surface, laser scattering by wind-generated foam as well as sun 
glint and solar scattering from the foam. An important feature of the algorithm is the use of the infrared 
backscattering as a basis for predicting the specular green backscattering. This allows the implicit 
correction for water surface geometry that may not follow simplistic model assumptions. 

Sensitivity analysis indicates that the model is insensitive to the wind speed estimate, but is 
relatively sensitive to errors in the atmospheric transmittance and total integrated attenuated 
backscatter. However, errors in the transmittance or total attenuated integrated backscatter that are 
positively correlated between the two bands will tend to cancel out because of the reversed sign of the 
sensitivity. 

In a case study, the method was applied to nighttime CALIOP data over Tampa Bay, using 
interpolated wind speed and visibility data from airports and buoys, and comparison was made to 
MODIS 645 mn remote sensing reflectance. The model was found to be very insensitive to the wind 
speed and visibility inputs. The results show a small but statistically significant correlation, and 
individual dates (particularly those with relatively high atmospheric transmittance) show much higher 
correlations. 

The CALIOP total integrated attenuated backscatter contains infonnation about scattering from the 
subsurface. More work may be needed to refine the models in order to obtain the highest signal-to- 
noise ratio possible from this algorithm. CALIOP was not designed to measure optical properties of 
the water column, but there is useful information hidden in these data. 
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